Nonequilibrium phase transition in a non integrable 
zero-range process 



C. Godrecheft 

f Isaac Newton Institute for Mathematical Sciences, 20 Clarkson Road, Cambridge, 
CB3 OEH, U.K. 

| Service de Physique de l'Etat Condense, CEA Saclay, 91191 Gif-sur-Yvette cedex, 
France 

Abstract. The present work is an endeavour to determine analytically features of 
the stationary measure of a non-integrable zero-range process, and to investigate the 
possible existence of phase transitions for such a nonequilibrium model. The rates 
defining the model do not satisfy the constraints necessary for the stationary measure 
to be a product measure. Even in the absence of a drive, detailed balance with respect 
to this measure is violated. Analytical and numerical investigations on the complete 
graph demonstrate the existence of a first-order phase transition between a fluid phase 
and a condensed phase, where a single site has macroscopic occupation. The transition 
is sudden from an imbalanced fluid where both species have densities larger than the 
critical density, to a critical neutral fluid and an imbalanced condensate. 
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1. Introduction 

Zero-range processes (ZRP) are minimal models [1], often used as simplified realisations 
of more complex processes (for reviews, see [2, 3, 4]). For instance they are instrumental 
for the understanding of condensation transitions in driven diffusive systems [5, 6]. They 
are closely related to urn models, which themselves are simplified models of a number 
of stochastic processes in Physics [8]. In a ZRP, particles, all equivalent, hop from sites 
to sites on a lattice, with prescribed rates which only depend on the occupation of the 
departure site. The fundamental property of this process is that the stationary measure 
is explicitly known as a function of the rates, and is a product measure [1, 7]. 

A natural generalisation of this definition is when two different species are allowed 
to coexist on each site, again hopping with prescribed rules. However in this case the 
stationary measure is a product measure only if the rates with which the particles of 
both species leave a given site satisfy a constraint [9, 10] (see eq. (2.5) below). For 
short, we refer to ZRP satisfying (2.5) as integrable. If these rates do not satisfy the 
constraints, the stationary measure is not known, and the corresponding ZRP is a generic 
nonequilibrium model: it violates detailed balance, even in the absence of a drive applied 
to the system, as will be shown below. A question of fundamental importance posed 
by the study of nonequilibrium systems is the nature of their stationary state, and in 
particular the possible existence of phase transitions at stationarity. 

The present work is devoted to the investigation of this question on the particular 
example of a non integrable two-species ZRP. The model arose from the study of a 
two-species driven diffusive system (DDS) exhibiting, at stationarity, condensation with 
coexistence between a high and a low density phase in each individual domain [6]. A 
domain in the original DDS, i.e. a stretch of particles of the two species, corresponds to 
a site in the ZRP, while the high and low density phases correspond to the two species 
of the ZRP. 

We study the model on the complete graph (i.e., in the fully connected geometry), 
using analytical and numerical methods. While for equal densities of the two species 
the transition between a fluid phase and a condensed phase is continuous (as is the 
case for the corresponding single-species ZRP), for non-equal densities this transition is 
discontinuous. The model exhibits a sudden phase transition from an imbalanced fluid 
where both species have densities larger than the critical density to a neutral fluid, with 
densities of both species equal to the critical density, and an imbalanced condensate. 
As a consequence reentrance is observed. The system is successively fluid, condensed, 
fluid, when increasing the density of one species, holding the density of the other species 
fixed. Coexistence between the two phases takes place along the transition line only. 
This study can serve as a template for the study of the one-dimensional model. 
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2. Definition of the model 

2.1. A reminder on zero-range processes 

We first give a short reminder of the definition of a ZRP. Consider, in any dimension, 
a lattice of M sites on which N particles are moving. Multiple occupancy of a site is 
allowed. The dynamics consists in choosing a site at random, then transferring one of 
the particles present on this site, to an arrival site. On the complete graph all sites are 
connected. The arrival site is any site chosen randomly. In one dimension, the arrival 
site is one of the two nearest neighbours, chosen with a given probability, p, to the right, 
or q = 1 — p, to the left. The transfer of particles is done with the rate u k (k > 0), only 
depending on the number k of particles on the departure site. 

The fundamental property of the ZRP is that its stationary measure is known, and 
is a product measure, as follows. Let us denote by iVj the random occupation of site %. 
The stationary weight of a configuration of the system is 

i M 

V(N 1 ,...,N M ) = - l[p Ni , (2.1) 

£m,n i= i 

where the normalisation factor Z M;N reads 

Zm,n = E • • • E m • • -Pnm 8 fe N i> N ) ■ (2-2) 

iVi N M V i / 

For a given rate u k , the factors pu obey the relation 

u k p k =p k _ 1 (2.3) 
which leads to the explicit form 

p = 1, p k = . (2.4) 

Ui...U k 

Let us emphasize two important characteristics of the ZRP (the same holding for 
integrable two-species ZRP, defined below). 

• When the dynamics is symmetric, e.g. in the one dimensional geometry with 
p = 1/2, or in the fully connected geometry, detailed balance with respect to 
the stationary measure is satisfied. For the single-species ZRP the detailed balance 
condition reads 

PkPlU k = p k -!Pi + iUi +1 , 

which is precisely the property that leads to (2.3)4 

• As can be seen on (2.1), the stationary measure is independent of the asymmetry. 
As a consequence, any property of the ZRP based on the sole knowledge of this 
measure is itself independent of the asymmetry. For example, with special choices 
of the transfer rate u k , a condensation transition can occur in the system. The 
features characterising this phase transition are independent of the asymmetry. 

| When the dynamics is not symmetric, the condition of detailed balance is replaced by a condition of 
pairwise balance. See [4] for a discussion of this point. 
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This is in contrast with the ZRP studied in the present work, where detailed balance 
is not satisfied, even when the dynamics is symmetric, i.e., in the absence of a drive, as 
explained below. In this sense this model is a generic nonequilibrium model, and the 
phase transition described in the next sections is specific of a nonequilibrium system. 

2.2. The model considered in the present work 

The model considered in the present work is a two-species ZRP. The general definition 
of a two-species ZRP is a simple extension of that of the usual ZRP [9, 10]. Consider, 
in any dimension, a lattice of M sites with n particles of type 1, m particles of type 

2. The dynamics consists in choosing a site at random, then transferring one of the 
particles present on this site, of one of the species chosen at random, to an arrival site. 
The transfer of particles is done with rates u k j (k > 0) for a particle of the first species, 
and v k j (/ > 0) for a particle of the other species, where k and I are respectively the 
number of particles of each species on the departure site. 

At variance with the case of single-species ZRP where the stationary measure is a 
product measure for any choice of the transfer rate, for a two-species ZRP this property 
holds only if the following constraint on the rates u k j and v k j is satisfied [9, 10]: 

u k ,iv k -i,i = Ufc,j«fc,i-i- (2.5) 

In the present work we choose rates which violates this constraint. As a consequence, 
nothing a priori is known on the nature of the stationary measure of the model. The 
rates read 

U k ,l = 1 + y, V k j = 1 + p (2.6) 

where b is a given parameter, which plays the role of inverse temperature [4]. We also 
set u k fi = v 0: i = 1 + b in order to complete the definition of the process. These rates 
favour equality of the number of particles of both species on each site. This choice aims 
at reproducing a feature of the original DDS, where inside the domains coexistence 
between a high and a low density phase takes place. The model was first introduced 
in [6], and studied in the equal density case. This is reviewed and extended below. We 
then focus on the non-equal density case. 

3. The case of two sites 

We begin by considering the case where the system is made of two sites. This case shares 
many common features with the complete model and serves as a useful preparation for 
the rest. 

A configuration of the system is entirely specified by the numbers k and I of particles 
of each species on site 1, since the number of particles on site 2 are then just equal to 
n — k and m — l. Therefore the weight of a configuration of the system is given by the 
probability f k ,i(t) that site 1 contains k particles of one species, and / particles of the 
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other species, at time t. It obeys the master equation 

^TT ' — u k+l,l — Sk,n) + ^fc.Z+l /fc,/+l(l — ^l.m) 

at 

+ U n _k+l,m-l fk-l,lO- ~ fikfi) + ^ra-fe.m-i+l /fc,J-l(l — <^,o) 

— [ M fc,Z + ffc,J + Un-k,m-l + fn-fc,m-j] /fc,Z, (3-1) 

where it is understood that w ,i = ffc,o — 0. This is the master equation of a biased 
random walk in the rectangle 0<k<n,0<l<m, with reflecting boundary conditions. 

We would like to know the stationary solution of this equation, for any 
choice of rates Uk,i, Vkj- It turns out that this question is already too hard to 
answer for the seemingly simple problem of a two-site model. We must content 
ourselves of the knowledge of the stationary solution for the class of processes fulfilling 
the constraint (2.5). Indeed, proviso the rates fulfill this constraint, the stationary 
distribution is given by 

e Pk,l Pn—k,m—l / Q \ 

Jk,l - 7y > \ 6 - 2 ') 



J M,n,m 



where 



Pk,lUk,l = Pk-l,l 

Pk,iVk,i = Pk,i-i, (3-3) 

and ZM, n ,m is a normalisation (the partition function). Relations (3.3) can be iterated, 
thus determining the pk,i in terms of the rates. Eqs. (3.2) and (3.3) generalize eqs. (2.1) 
and (2.3). 

The method used in [9, 10] to obtain these results consists in making the 
ansatz (3.2), carry this form into the master equation, which leads to (3.3), which 
itself imposes (3.5) compatibility relation. 

We wish to bring an independent and complementary viewpoint to this issue. We 
first note that the dynamics between the two sites is symmetric. We therefore question 
the possibility for the process to be reversible in time, and the consequences thereby. 
Reversibility is equivalently the property that the process obeys detailed balance with 
respect to the stationary measure, or otherwise stated that the system is at equilibrium. 
We proceed as follows. 

(i) We first determine the stationary distribution when detailed balance is obeyed. 
Consider the transitions from {k, 1} to {k + 1, 1} and back, and from {k, 1} to {k, I + 1} 
and back. Detailed balance requires 

U n -k,m~l fk,l — Uk+l t l fk+l,l-i 
V n -k,m-l fk,l = Vk,l+1 fk,l+l- 

It is readily found that a solution of these equations is given by (3.2) and (3.3). 

(ii) We now determine, by yet another path, the conditions on the rates for the model to 
satisfy reversibility. We use the Kolmogorov criterion [11, 12] which is a necessary and 
sufficient condition for a Markov process to be reversible. This condition states that the 
product of the transition rates along any cycle in the space of configurations should be 
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equal to the product of the transition rates for the reverse cycle. In the present case, 
the space of configurations is the rectangle < k < n, < I < m. Taking the cycle 

(k, I) -> (k, I - 1) -> (k + 1, 1 - 1) -> (k + 1, /) -> (A;, /), 

then the cycle in reverse order, the Kolmogorov condition leads to the equation 

Uk+l,lVk,l _ U n -k,m-l Vn-k,m-(l-l) ^ ^ 

U n -k,m-(!-i)V n -(k+i),m-(l-l) 

The two sides of this equation should be satisfied independently. This imposes that 

Uk,i Vk-i,i = u k ,i-i v k ,i, (3.5) 

which is the constraint (2.5) .§ 

To summarize, reversibility implies stationary product measure, eqs. (3.2), (3.3), 
and a constraint on the rates, eq. (3.5). The reciprocal statement holds. The proof 
follows easily from the fact that a Markov process with a finite configuration space has 
a unique stationary solution. We leave it to the reader. 

The physical interpretation of the results above is that when the system is at 
equilibrium, its energy is equal to the sum of the energies of two independent sites. 

Conversely, for a choice of rates violating (3.5), as is the case for the model studied 
here, the model is not reversible, the stationary measure does not take the simple form 
(3.2), and is not known a priori. In other words, for a general choice of rates, the two-site 
model can have an arbitrarily complex stationary measure. In this sense it represents 
an example of a minimal nonequilibrium system. 



4. The model on the complete graph 

The virtue of considering the fully connected geometry, in the thermodynamical limit 
of an infinite system, is that it leads to analytical results on the model. The rest of the 
paper is devoted to this case. 

Consider again the single-site occupation probability fk,i(t), that is the probability 
that a generic site contains k particles of one species, and / particles of the other species, 
at time t. Conservation of probability and particle numbers imposes J2k,i fk,i(t) = 1 and 

oo oo 

E */*(*) = ah, EW) = pa, 

k=l 1=1 

where for a large system, densities are defined as p 1 = n/M, p 2 = m/M, and where the 

marginals are denoted by f k = £; f k ,i, and f x = £ fc fk,i- 

§ Eq. (3.4) can be satisfied by imposing symmetry relations on the rates: 

Uk+l,l = U n -k, m -l,Vk,l+l = V n -k, m -l- 

The corresponding stationary measure is uniform: 

f 1 

m (n + l)(m+l)' 

and detailed balance is obeyed. We discard this solution because the rates would then also depend on 
the arrival site. 
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The master equation for the temporal evolution of fk,i{t) reads 
dfkM 



Uk+l,l fk+l,l + Vk,l+l fk,l+l 

u t /fc-i,z(l - 4,o) + v t 
[Uk,l + v k ,i + u t + V t ] fk,< 



(it 

+ u t / fc _i,j(l - 4,o) + vt fk,i-i(l - &i,o) (4-2) 



where 



U t = U k,l fk,l, v t = J2 V k,l fk,l 

k . I k .1 



are the mean rates at which a particle arrives on a site (k — > k + 1) or (I — > I + 1). 
Equation (4.2) is the master equation for a biased random walk in the quadrant k, I > 0, 
with reflecting boundary conditions on the axes. 

We wish to determine the stationary solution of this equation. We follow the same 
line of thought as in the previous section. We show that the stationary distribution f kji 
has a known closed-form expression only if reversibility is assumed. Indeed, using the 
detailed balance conditions 

fk+1,1 u k+i,i — u fk,l, 
fk,l+l v k,l+l — v fk,h 

it is easy to derive the following explicit expression for the stationary distribution f k y. 
_ p Kl u k v l 

Jk,i — ^ =nr=7' 

Y,k,lPk,lU k V 1 

where the pkj are given by (3.3), and u and v are the stationary mean hopping rate. 

Let us also show that, as for the two-site system, the constraint (3.5) is a 
consequence of imposing reversibility. Indeed, the space of configurations is the quadrant 
k, I > 0. Taking the cycle 

(k, I) -> (h, I - 1) -> (k + 1, 1 - 1) -> (k + 1, 1) -> (k, /), 

then the cycle in reverse order, the Kolmogorov condition implies 

V kt iUVU k+1 j = UVk+^iUk+xj-xV, 

which yields (3.5). 

As for the two-site system, we conclude that conversely, when (3.5) is violated, as 
is the case with the choice of rates (2.6), the stationary distribution remains unknown. 
The present work is an endeavour to determine features of the stationary measure of 
the model for a thermodynamical system, and to investigate the possible existence of 
nonequilibrium phase transitions. 

In the case of an integrable two-species ZRP, the fugacities are functions of the 
densities (see eq. (4.3)). The duality fugacity-density is replaced here by the duality 
mean hopping rate-density. 
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5. Criticality 

As will appear clearly as we proceed, the critical point for this model is unique, and 
corresponds to taking u — v — 1. 



5.1. Continuum limit: universal properties 

Let us first consider the continuum limit of the stationary equation, in the asymptotic 
regime where k and I are large, and setting u = e M , v = e u , where /i and v are small. 
Expanding fki to second order, we obtain 



dk 2 



dl 2 



dk \l 



dl \k 



(5-1) 
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Figure 1. Decay exponent a as a function of b. At large values of b, a ~ 26+\/2. 



At criticality, i.e. when /i = z/ = 0, eq. (5.1) becomes scale invariant, and reads 
d 2 fk,i d 2 f k i (ldfki ldf k i\ 

-cW + ^ +b {i^k- + k^r) =0 - (5 - 2) 

Using polar coordinates: k = rcosO, I = rsin^, with < 9 < n/2, this equation is 
transformed into 

^ + 4^ + I^f 1 + ^).0. (5.3, 
or 2 r l oo z r or \ sin 26 J 

Now, setting f(r,6) = r~ a g(8), we find an equation for the angular function g(9): 

The unknown decay exponent a is determined by the boundary conditions imposed on 
g(0), which are the quantisation conditions of this Schrodinger equation. Indeed, g(9) 
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is positive for < 6 < ir/2, symmetric with respect to 7r/4 and must vanish for 9 = 
or 6 = 7r/2. For special values of b, exact solutions of eqs. (5.3) can be found: 

f(r,6) = r- 2 sin 6 cos 6 (b = 0), (5.5) 

f(r,6) = r- 3 sin6 cos6(sm6 + cos6) (5 = 2/3). (5.6) 

For b = the original model has no critical behaviour, hence formally a = 0. On the 
other hand the prediction of the continuum limit for the decay exponent, in the limit 
b — > 0, is a = 2. The decay exponent a is discontinuous at b = 0. 

For a generic value of b, the decay exponent a is determined by numerical integration 
of the differential equation (5.4). At large values of b, the behaviour of this exponent 
can be obtained analytically. Indeed, expanding the potential term in (5.4) to second 
order around its minimum, located at 7r/4, yields the equation of a harmonic oscillator, 
with coupling constant oj = 2\fah and energy a(a — 2b) /2: 

g"{x) + a(a -2b- 4bx 2 )g{x) = 

where x = 7r/4 — 6. Imposing that the ground state energy be equal to u/2 yields the 
asymptotic quantisation condition a = 2b + \[2. (See figure 1.) 

As a consequence, we find the remarkable result that, at criticality, the marginal 
distributions fk and fi decay as power laws at large occupations, with a non-trivial 
exponent equal to a — 1. The same holds for p m , with to = k + 1, the distribution of the 
total number of particles on a site, 

m 

Pm = Yl fk,m-k ~ m'^"^. 
k=0 

Note finally that both the function g(8) and the exponent a are universal, and only 
depend on b. 

5.2. Discrete equations: critical density 

The determination of the non universal critical density p c of both species, where 

oo oo 

Pc= £fc/k = (U = V = 1), 

k=l 1=1 

requires the knowledge of the stationary solution of the discrete eqs. (4.2). These 
are integrated numerically with u — v — 1, using the following method. We truncate 
these equations at a given value of k + I denoted by to*, which plays the role of a 
cut-off. We solve the linear system AF = I, where F is the column matrix of the 
occupation probabilities fkj, I is the matrix containing the inhomogeneous term / 00 , 
itself determined at the end of the computation by normalisation, and A is the matrix 
deduced from the stationary equations. We impose the boundary conditions = 
outside the triangle delimited by k — 0, I — 0, and k + 1 — to*. The maximal value of 
the cut-off to* attainable is limited by the size of the matrices involved. For example, 
taking m* = 160 corresponds to a linear system of order 13040. 
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As an illustration we take b = 3/2, corresponding to a value of the decay exponent 
a « 4.520. Extrapolating the data for several values of m*, using the estimate 

p c - p c (m+) m / dmmp m ~ m~ (a_3) , 



as depicted in figure 2, leads to p c ~ 0.976. 
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Figure 2. Determination of the critical density by extrapolation of the data for 
m* = 40, 60, . . . , 160. The circle on the vertical axis is the extrapolated value 
for p c . (b = 3/2, a = 4.520 ...,u = v = l.) 



The theoretical prediction for the critical decay exponent of p m , or f k = f h agrees 
perfectly well with numerical measurements. 

6. Fluid phase 

For non zero values of p and v the system is driven away from criticality. We begin 
by investigating exponentially decreasing solutions of the continuum limit stationary 
equation (5.1). We then study those of the discrete stationary equations. We finally 
determine the region of existence of such solutions. 

6.1. Stationary solutions at exponential order: continuum limit 

If we content ourselves of the knowledge of the stationary solutions at exponential order, 
that is retaining their exponential dependence only, and discarding any prefactor, then 
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terms containing b can be neglected. Equation (5.1) now reads 

d 2 fk,i . d 2 fk,i dfk,i dfk,i , s 

Setting / M = e ^ k+ ^ h k ,i in (6.1) yields 

d 2 h K i d 2 h kyl /i 2 + ^ 2 , _ 

~cW~ ~~dP — hh * ~ °" 

Changing to polar coordinates and setting h(r, 9) = u(r)v(9), we obtain, after rescaling 
r by vV 2 + z/2 /2, 

rV(r) + ru'(r) ~ ( r2 + n 2 )u(r) = 0, (6.2) 
where n is to be determined, and 
v"(9)+n 2 v(9) = 0. 

Imposing v(9) = for 9 = and 7r/2 leads to i>(0) = sin 29, hence n — 2. The solution 
of the differential equation for w(r) is the Bessel function 

u(r) = K 2 (V/i 2 + v 2 • 

Finally, the solution, when 6^0, reads, up to a normalising constant, 

r 



f(r,9) = K 2 ^/i 2 + z/ 2 -J e f(Mcose+,sine) gin2 ^ ( 6 ^ 

This solution encompasses all three regimes where /ir ~ 1, /ir <C 1, or /ir >> 1. Simplified 
expressions are obtained in the two latter cases: 

• For small values of the argument, we have K 2 (x) ~ 1/x 2 . We thus find 

f(r, 9) const, r~ 2 sin 26* , (/ir < 1), 

which matches consistently with the critical solution (5.5) f or b — > 0. 

• For large values of the argument, we have K 2 (x) as ^Jir/2xe~ x , hence we obtain 
the asymptotic behaviour 

/(r, 0) « const. r~ 1/2 e~ rP(e) sin 20, (/ir > 1), (6.4) 

where 

1 



P(0) = - \^Jfj? + z/2 - fj, cos - v sin #J . (6.5) 

For any values of /i and v non simultaneously positive, P{9) is positive, / is exponentially 
decaying, corresponding to a fluid phase. When \x and z/ are simultaneously positive, 
P{9) vanishes at an angle 9 satisfying tan 6* = v/n- For such a value of 9, the function 
f(r, 9) ~ r -1 / 2 is not normalisable. The whole region /i > and z/ > is therefore non 
physical. 
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6.2. Stationary solutions at exponential order: discrete equations 

In order to investigate exponentially decaying solutions beyond the continuum limit, we 
consider again the discrete stationary equations. As above, for k and I large, we can 
neglect terms containing b, thus obtaining 

fk+1,1 + fk,i+i + u fk- lt i + v - (2 + u + v)f kjl = 0. (6.6) 

Introducing the generating function f{x,y) = fk,ix k y l , we get from (6.6) 

D(x,y)f(x,y) = A(x,y), (6.7) 

where 

D(x, y) = x- 1 +y- 1 -2 + u(x - 1) + v(y - 1). (6.8) 

The locus of singularities of / is thus given by D(x, y) = 0. The right-hand side, A(x, y), 
comes from the contribution of the boundary terms f(0,y) and f(x, 0) ||. By inversion 
we have 

h* = ip-p-fav)*- k v- 1 - 

J 2m x 2my 

At large k and /, can be estimated by taking the saddle point of this expression, 
yielding 

where x, y is a point on the curve D(x,y) = 0, and P{0) = cos 6 In x + sin^lny. 
Extremising P{6) on this curve with respect to the variables x and y, i.e. the expression 
P{6) — \D(x,y), where A is a Lagrange multiplier, leads, together with D(x,y) = 0, to 
three equations, which determine A, x and y, hence P(0). 

Let us first check that this method leads to the expected result (6.5) in the particular 
case of the continuum limit. Set x = e s , y = e*, u = e M , v = e u , with fi and v small. 
The equation for D(x,y) reads: 

s 2 + t 2 + lis + vt = 0. 

Extremising P{6) — \D(x,y) with respect to s and t yields 

A cos# — 2s + /i = 0, 
\sm9-2t + v =0. 

The three former equations lead, after some algebra, to eq. (6.5). 

The general case leads to lengthy expressions for P(9). We give the results of this 
method for a particular example. We choose u = 1.89, v = 0.66 in the stationary 
equations, which we integrate by solving the linear system, as explained in section 5.2, 

|| For the generic case b > 0, A(x,y) is finite on D(x,y), while if b = 0, the singularities of A(x,y) 
cancels those of D(x,y) in such a way that the resulting expression for ,f(x,y) has just a simple pole 
in both complex variables x and y: 

(1 - ux)(l - vy) 
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for increasing values of the cut-off m*. The resulting values of the densities p\ and p 2 
are plotted in figure 3. Clearly p 2 is larger than p c V. Figure 4 depicts hxp m , as obtained 
by the same method, together with the theoretical prediction for the coefficient of the 
exponential decay: 



Pr, 



J rdrd9e- rP{0) 5 



m 



cos 8q +sm f 



cos 9 + sin 9 , 

where 9q denotes the value of the angle such that the argument of the exponential is 
minimum. In the present case, 9q = 0, and p m ~ e -o.o372m_ 
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Figure 3. Densities as functions of m+ (u = 1.89, v = 0.66, b = 3/2). The 
dashed line corresponds to p c . 



6. 3. Domain of existence of the fluid phase in the u — v plane 

The domain of existence of the homogeneous fluid solution in the u — v plane is shown in 
figure 5. It is the interior of the domain delimited by the two symmetric curves. These 
curves are obtained as follows. Consider the situation where one of the densities, p\, 
say, is infinite. Then v^j is to be taken equal to 1, and it is intuitively expected that 
the two species decouple. Hence fi — (1 — v)v l . It follows that 

v 

1 — V 

and 

u = Y,fi[l + 1 )=l + b(l-v)(l- ln(l - v)). (6.9) 



The two former equations give the equation of the boundary of the domain of existence 
of fluid solutions, for p 2 varying. The second part of the curve is obtained symmetrically 
by doing the same analysis with p 2 infinite. We note that the tangents to the two curves 
at the symmetric point u — v — 1 are parallel to the axes. 

V The values of u and v were precisely chosen to serve this purpose. We first integrated the master 
equation (4.2) numerically for b = 3/2, corresponding to p c ~ 0.976, and for values of the densities 
pi = 10 and p2 = 1. The stationarity values u « 1.894, v « 0.661 were thus obtained. 
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Figure 4. Distribution of the total single-site occupation p m as a function of 
m (m < m± = 140). Line: theoretical prediction for the coefficient of the 
exponential decay, (u = 1.89, v = 0.66, b = 3/2.) 



6.4- Limits of stability of the fluid phase 

Finally the question is how the domain of existence of the fluid region defined above 
is mapped onto the density plane. As we now show, this domain maps onto a region 
of the density plane complementary to a wedge with tip located at the critical point 
P\= Pi = Pci as depicted in figure 6. All the analysis relies on how the neighbourhood 
of the point u — v — 1 is mapped onto the density plane. 

We begin by a linear analysis of the mapping between the critical points (u = v = 1) 
and (pi = p2 = p c )- Consider a small segment, with one of the two ends located at 
u — v — 1, and with the other one at a given angle with the u axis. Let t be the 
tangent of this angle. Because of the symmetry between the two species, and since the 
transformation of the local derivatives around the critical point is linear, the slope of 
the transformed segment in the density plane is given by 

T=l±^. (6.10) 
c + t 

Thus, if t = ±1, then T = ±1. The constant c is determined numerically by taking 
t = oo (segment parallel to the v axis). The limiting slopes of the tangents to the wedge 
at the tip follow from (6.10). For the lower edge it is given by T(t = 0), i.e., 1/c. For 
6 = 3/2, we find T(t = 0) w 0.48, with m* = 80. 

We then investigate how points located on the boundary curve (6.9), at increasing 
distances of the critical point, u = v = 1, are mapped onto the density plane. For 
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Figure 5. Domain of existence of the homogeneous fluid solution (b = 3/2). 
The lower right wing corresponds to p\ = do, P2 finite, and symmetrically for 
the left upper wing. 

successive values of the cut-off m* the images of these points are expected to converge 
to a single curve. The lower edge of the wedge depicted in figure 6 is the curve with 
m* = 160. The other edge is obtained by symmetry. These edges represent the limits 
of stability of the fluid phase. 

7. Condensation and phase diagram 

7.1. Equal densities 

The case of equal densities, p\ — p%, is similar to the situation encountered for a single- 
species ZRP [3, 4]. From the analysis of the previous sections, as well as from numerical 
integrations of the temporal equations (4.2), or of their stationary form, the following 
picture is obtained. 

The region u = v < 1 maps onto the fluid phase p\ = p 2 < p c , corresponding 
to exponential solutions of eq. (5.1). The critical point corresponds to u — v — 1, 
i.e., pi = p 2 = Pc- Condensation occurs for a > 3, i.e. b > 2/3 (see eq. (5.6)), and 
Pi = P2 > Pc- A condensate appears sustaining the excess density with respect to the 
critical fluid. The region u = v > 1 is unphysical. 
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7.2. Nonequal densities: Existence of a line of transition 

The limits of stability of the condensed phase are given by the line p 2 = p c , Pi > Pc 
and the symmetric line with respect to the bisectrix. 

There is numerical evidence for the existence of a transition line between the fluid 
phase and the condensed phase, lying in between the two corresponding stability lines 
of these phases. The transition is discontinuous on this coexistence line. There is no 
coexisting solutions accessible dynamically on both sides of the line. 

The transition between the fluid and condensed phases is obtained by Monte Carlo 
simulations of the model, using the following procedure. The density p\ is fixed to a 
given value greater than p c , and p 2 increases from a value less than p c . Crossing the 
stability line of the condensed phase, i.e. for p 2 > p c , one might expect condensation to 
occur. Instead, the only accessible phase turns out to be the fluid one (see an example 
of fluid solution in section 6.2). Then, increasing the density p 2 , and crossing the 
transition line, there is a sudden phase transition from an imbalanced fluid where both 
species have densities larger than the critical density, to a neutral critical fluid and an 
imbalanced condensate. Beyond this line, the only accessible solutions are condensed, 
with u — v — 1, while fluid solutions to eq. (4.2) do exist, as long as p 2 has not reached 
the edge of the wedge. A surprising consequence of this phase diagram is the occurrence 
of a reentrance phenomenon: increasing p 2 beyond the symmetric transition line (with 
respect to the bisectrix) the system becomes fluid again. 

We now describe more precisely the method for the determination of the location 
of the transition line. We fix the value of pi, for example p\ = 10, and let p 2 increase 
from a value less than p c . Then u and v are measured. Focussing on u, we observe that, 
when p 2 crosses some value, p 2 ~ 1.8, there is a sudden discontinuity in u, dropping from 
u ~ 1.4 down to u ~ 1. More precisely, for a system of given size, say M = 40, 60, . . ., 
the system is first run to stationarity. Then n successive runs of duration At less 
than the flipping time r between the fluid and the condensed phases, and such that 
nAt 3> t are performed. This flipping time is measured to be exponentially increasing 
with the system size (it is approximately doubled when M is incremented by 20). The 
histogram of the values of u is a bimodal distribution. The criterion for the location of 
the transition point (pl,p 2 ), when p 2 varies, consists in choosing the value of p 2 such 
that the two weights of the two maxima are equal. Thus for M = 40, 60, 80, we have 
p 2 ~ 1.86, 1.82, 1.8, respectively. 

We proceed in the same fashion to obtain the transition points visible on figure 6, 
for pi = 4, 6, 8. As p\ decreases down to p c , the discontinuity in u is smaller and the 
determination of the transition point is harder since it involves larger system sizes. 

8. Final remarks 

Let us first summarize the main outcomes of the present work. For the process 
considered, a non-integrable two-species ZRP, we are able to obtain a number of 
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Figure 6. Phase diagram in the density plane. Dot-dashed line with circles: 
line of transition points (the symmetric line is not figured). Dashed lines: limits 
of stability of the fluid phase. Dotted lines: limits of stability of the condensed 
phase. Straight lines at the tip (critical point) are the local tangents, computed 
as explained in the text, (b = 3/2, p c 0.976.) 

analytical results from the study of the model on the complete graph. In particular 
the critical phase is well understood. The coefficient of the exponential decay of the 
fluid solutions is analytically predicted. Finally we can predict the phase diagram of 
the system by a joint analytical and numerical investigation. A salient feature of the 
phase diagram is the presence of a single critical point. A more thorough analysis of 
the first-order phase transition taking place between the fluid phase and the condensed 
phase would be interesting, though probably hard to achieve. The existence of such a 
transition is not expected in integrable two-species ZRP [9, 10, 13, 14]. 

Beyond the present work, a natural question to ask is whether the phenomena 
observed for the fully connected geometry survive in the one-dimensional geometry. 
Preliminary investigations using the analysis performed on the complete graph as a 
template indicate similar behaviour. 
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